www.gusucode.com > 循环自相关函数工具箱源码程序 > matlab代做 修改 程序循环自相关函数工具箱/cyclostationary_toolbox/Contents.m
cyclostationary_toolbox/cyclic_3rd_order_cumulant.m 0100644 0031654 0031654 00000003034 06436565164 0023325 0 ustar 00acmc acmc 0000304 0001726 function C3=cyclic_3rd_order_cumulant(x1,x2,x3,alpha,max_tau)
%
% CYCLIC_3RD_ORDER_CUMULANT
%
% calculates the cyclic third order cumulant of
% three signals x1,x2,x3 at frequency alpha
%
% C3(k*alpha,tau1,tau2)=E{(x1(t)-E{x1(t)}) *
% (x2(t+tau1)-E{x2(t+tau1)} *
% (x3(t+tau2)-E{x3(t+tau2)} *
% exp(-jk(alpha)t) }
% for k=0 ... 1/alpha
%
%
% USAGE
% C3=cyclic_3rd_order_cumulant(x,y,alpha,max_tau)
%
% File: cyclic_3rd_order_cumulant.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=5
error('Incorrect number of arguments for function cyclic_3rd_order_cumulant');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2 pi in function cyclic_3rd_order_cumulant');
end
% Remove cyclic mean from signals
cmx1=cyclic_mean(x1,alpha);
cmx2=cyclic_mean(x2,alpha);
cmx3=cyclic_mean(x3,alpha);
lx=length(x1);
t=0:lx-1;
T=ceil(2*pi/alpha)-1;
for k=1:lx
x1(k)=x1(k)-1/(2*pi)*sum(cmx1.*exp(j*alpha*(0:T)*(k-1)));
x2(k)=x2(k)-1/(2*pi)*sum(cmx2.*exp(j*alpha*(0:T)*(k-1)));
x3(k)=x3(k)-1/(2*pi)*sum(cmx3.*exp(j*alpha*(0:T)*(k-1)));
end
C3=zeros(max_tau,max_tau,T+1);
ix=1:lx-max_tau-1;
for tau1=0:max_tau
for tau2=0:max_tau
for k=0:T
C3(tau1+1,tau2+1,k+1)=mean(x1(ix).*x2(tau1+ix) ...
.*x3(tau2+ix).*exp(j*k*alpha*t(ix)));
end
end
end
cyclostationary_toolbox/cyclic_3rd_order_cumulant_fast.m 0100644 0031654 0031654 00000006074 06436555370 0024347 0 ustar 00acmc acmc 0000304 0001726 function C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,T,max_tau)
%
% CYCLIC_3RD_ORDER_CUMULANT_FAST
%
% calculates the cyclic third order cumulant of
% three signals x1,x2,x3 at frequency alpha using a fast
% approximation based on the synchronous average of the
% time varying autocorrelation. Fundamental signal
% period can be defined as a single period or as a sequence
% of once per period pulse times.
%
% C3(k*alpha,tau1,tau2)=E{(x1(t)-E{x1(t)}) *
% (x2(t+tau1)-E{x2(t+tau1)} *
% (x3(t+tau2)-E{x3(t+tau2)} *
% exp(-jk(alpha)t) }
% for k=0 ... 1/alpha
%
%
% USAGE
% C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,alpha,max_tau)
%
% File: cyclic_3rd_order_cumulant_fast.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=5
error('Incorrect number of arguments for function cyclic_third_order_cumulant_fast');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function cyclic_third_order_cumulant_fast');
end
[rows,cols]=size(x1);
if rows>cols
x1=x1';
end
[rows,cols]=size(x2);
if rows>cols
x2=x2';
end
[rows,cols]=size(x3);
if rows>cols
x3=x3';
end
% Calculate synchronous averages from signals
mx1=synchronous_average(x1,T);
mx2=synchronous_average(x2,T);
mx3=synchronous_average(x3,T);
% Remove excess samples due to non-integer sampling
% and renove cyclic mean from signal
if length(T)==1
cp=1;np=1;
while cp+T<length(x1)
x1(cp:cp+floor(T)-1)=x1(cp:cp+floor(T)-1)-mx1;
x2(cp:cp+floor(T)-1)=x2(cp:cp+floor(T)-1)-mx2;
x3(cp:cp+floor(T)-1)=x3(cp:cp+floor(T)-1)-mx3;
cp=cp+floor(T);
np=np+T;
if (np-cp)>1
x1=[x1(1:cp-1);x1(cp+1:length(x1))];
x2=[x2(1:cp-1);x2(cp+1:length(x2))];
x3=[x3(1:cp-1);x3(cp+1:length(x3))];
np=np-1;
end
end
n=floor((length(x1)-2*max_tau-1)/T);
else
% extract time series correlated with periodic pulses
rot_positions=T;
T=floor(median(diff(rot_positions)));
nx1=[];
nx2=[];
nx3=[];
n=length(rot_positions)-2;
for k=1:n;
cp=rot_positions(k);
nx1=[nx1; x1(cp:cp+T-1)-mx1];
nx2=[nx2; x2(cp:cp+T-1)-mx2];
nx3=[nx3; x3(cp:cp+T-1)-mx3];
end
nx1=[nx1;x1(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx1(1:tau+1)];
x1=nx1;
nx2=[nx2;x2(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx2(1:tau+1)];
x2=nx2;
nx3=[nx3;x3(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx3(1:tau+1)];
x3=nx3;
end
% Compute time varying third order cumulant
r=zeros(max_tau+1,max_tau+1,floor(T));
temp=zeros(floor(T),n);
t=(1:floor(T)*n);
for tau1=0:max_tau
for tau2=0:max_tau
temp(:)=x1(t).*x2(t+tau1).*x3(t+tau2);
r(tau1+1,tau2+1,:)=mean(temp');
end
end
% Take DFT of time varying toc
C3=zeros(max_tau+1,max_tau+1,floor(T));
for tau1=0:max_tau
for tau2=0:max_tau
C3(tau1+1,tau2+1,:)=fft(r(tau1+1,tau2+1,:))/T;
end
end
cyclostationary_toolbox/cyclic_4th_order_cumulant.m 0100644 0031654 0031654 00000004173 06436565216 0023337 0 ustar 00acmc acmc 0000304 0001726 function C4=cyclic_4th_order_cumulant(x1,x2,x3,x4,alpha,max_tau)
%
% CYCLIC_4TH_ORDER_CUMULANT
%
% calculates the cyclic fourth order cumulant of
% three signals x1,x2,x3,x4 at frequency alpha
%
% C3(k*alpha,tau1,tau2,tau3)=E{(x1(t)-E{x1(t)}) *
% (x2(t+tau1)-E{x2(t+tau1)} *
% (x3(t+tau2)-E{x3(t+tau2)} *
% (x4(t+tau3)-E{x4(t+tau3)} *
% exp(-jk(alpha)t) }
% for k=0 ... 1/alpha
%
%
% USAGE
% C3=cyclic_4th_order_cumulant(x,y,alpha,max_tau)
%
% File: cyclic_4th_order_cumulant.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=6
error('Incorrect number of arguments for function cyclic_4th_order_cumulant');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2 pi in function cyclic_4th_order_cumulant');
end
% Remove cyclic mean from signals
cmx1=cyclic_mean(x1,alpha);
cmx2=cyclic_mean(x2,alpha);
cmx3=cyclic_mean(x3,alpha);
cmx4=cyclic_mean(x4,alpha);
lx=length(x1);
t=0:lx-1;
T=ceil(2*pi/alpha)-1;
for k=1:lx
x1(k)=x1(k)-1/(2*pi)*sum(cmx1.*exp(j*alpha*(0:T)*(k-1)));
x2(k)=x2(k)-1/(2*pi)*sum(cmx2.*exp(j*alpha*(0:T)*(k-1)));
x3(k)=x3(k)-1/(2*pi)*sum(cmx3.*exp(j*alpha*(0:T)*(k-1)));
x4(k)=x4(k)-1/(2*pi)*sum(cmx4.*exp(j*alpha*(0:T)*(k-1)));
end
C4=zeros(max_tau,max_tau,max_tau,T+1);
ix=1:lx-max_tau-1;
for tau1=0:max_tau
for tau2=0:max_tau
for tau3=0:max_tau
for k=0:T
M4=mean(x1(ix).*x2(tau1+ix).*x3(tau2+ix) ...
.*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
c2_12=mean(x1(ix).*x2(tau1+ix).*exp(j*k*alpha*t(ix)));
c2_34=mean(x3(tau2+ix).*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
c2_13=mean(x1(ix).*x3(tau2+ix).*exp(j*k*alpha*t(ix)));
c2_24=mean(x2(tau1+ix).*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
c2_14=mean(x1(ix).*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
c2_23=mean(x2(tau1+ix).*x3(tau2+ix).*exp(j*k*alpha*t(ix)));
C4(tau1+1,tau2+1,tau3+1,k+1)=M4-c2_12*c2_34-c2_13*c2_24-c2_14*c2_23;
end
end
end
end
cyclostationary_toolbox/cyclic_4th_order_cumulant_fast.m 0100644 0031654 0031654 00000007617 06436562676 0024371 0 ustar 00acmc acmc 0000304 0001726 function C4=cyclic_4th_order_cumulant_fast(x1,x2,x3,x4,T,max_tau)
%
% CYCLIC_4TH_ORDER_CUMULANT_FAST
%
% calculates the cyclic fourth order cumulant of
% four signals x1,x2,x3,x4 at frequency alpha using a fast
% approximation based on the synchronous average of the
% time varying autocorrelation. Fundamental signal
% period can be defined as a single period or as a sequence
% of once per period pulse times.
%
% C4(k*alpha,tau1,tau2)=E{(x1(t)-E{x1(t)}) *
% (x2(t+tau1)-E{x2(t+tau1)} *
% (x3(t+tau2)-E{x3(t+tau2)} *
% (x4(t+tau2)-E{x4(t+tau2)} *
% exp(-jk(alpha)t) }
% for k=0 ... 1/alpha
%
%
% USAGE
% C4=cyclic_4th_order_cumulant_fast(x1,x2,x3,x4,alpha,max_tau)
%
% File: cyclic_4th_order_cumulant_fast.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=6
error('Incorrect number of arguments for function cyclic_4th_order_cumulant_fast');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function cyclic_4th_order_cumulant_fast');
end
[rows,cols]=size(x1);
if rows>cols
x1=x1';
end
[rows,cols]=size(x2);
if rows>cols
x2=x2';
end
[rows,cols]=size(x3);
if rows>cols
x3=x3';
end
[rows,cols]=size(x4);
if rows>cols
x3=x3';
end
% Calculate synchronous averages from signals
mx1=synchronous_average(x1,T);
mx2=synchronous_average(x2,T);
mx3=synchronous_average(x3,T);
mx4=synchronous_average(x4,T);
% Remove excess samples due to non-integer sampling
% and renove cyclic mean from signal
if length(T)==1
cp=1;np=1;
while cp+T<length(x1)
x1(cp:cp+floor(T)-1)=x1(cp:cp+floor(T)-1)-mx1;
x2(cp:cp+floor(T)-1)=x2(cp:cp+floor(T)-1)-mx2;
x3(cp:cp+floor(T)-1)=x3(cp:cp+floor(T)-1)-mx3;
x4(cp:cp+floor(T)-1)=x4(cp:cp+floor(T)-1)-mx4;
cp=cp+floor(T);
np=np+T;
if (np-cp)>1
x1=[x1(1:cp-1);x1(cp+1:length(x1))];
x2=[x2(1:cp-1);x2(cp+1:length(x2))];
x3=[x3(1:cp-1);x3(cp+1:length(x3))];
x4=[x4(1:cp-1);x4(cp+1:length(x4))];
np=np-1;
end
end
n=floor((length(x1)-2*max_tau-1)/T);
else
% extract time series correlated with periodic pulses
rot_positions=T;
T=floor(median(diff(rot_positions)));
nx1=[];
nx2=[];
nx3=[];
nx4=[];
n=length(rot_positions)-2;
for k=1:n;
cp=rot_positions(k);
nx1=[nx1; x1(cp:cp+T-1)-mx1];
nx2=[nx2; x2(cp:cp+T-1)-mx2];
nx3=[nx3; x3(cp:cp+T-1)-mx3];
nx4=[nx4; x4(cp:cp+T-1)-mx4];
end
nx1=[nx1;x1(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx1(1:tau+1)];
x1=nx1;
nx2=[nx2;x2(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx2(1:tau+1)];
x2=nx2;
nx3=[nx3;x3(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx3(1:tau+1)];
x3=nx3;
nx4=[nx4;x4(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx4(1:tau+1)];
x4=nx4;
end
% Compute time varying third order cumulant
r=zeros(max_tau+1,max_tau+1,max_tau+1,floor(T));
temp=zeros(floor(T),n);
t=(1:floor(T)*n);
for tau1=0:max_tau
for tau2=0:max_tau
for tau3=0:max_tau
temp(:)=x1(t).*x2(t+tau1).*x3(t+tau2).*x4(t+tau3);M4=mean(temp');
temp(:)=x1(t).*x2(t+tau1);c2_12=mean(temp');
temp(:)=x3(t+tau2).*x4(t+tau3);c2_34=mean(temp');
temp(:)=x1(t).*x3(t+tau2);c2_13=mean(temp');
temp(:)=x2(t+tau1).*x4(t+tau3);c2_24=mean(temp');
temp(:)=x1(t).*x4(t+tau3);c2_14=mean(temp');
temp(:)=x2(t+tau1).*x3(t+tau2);c2_23=mean(temp');
r(tau1+1,tau2+1,tau3+1,:)=M4-c2_12.*c2_34-c2_13.*c2_24-c2_14.*c2_23;
end
end
end
% Take DFT of time varying foc
C4=zeros(max_tau+1,max_tau+1,max_tau,floor(T));
for tau1=0:max_tau
for tau2=0:max_tau
for tau3=0:max_tau
C4(tau1+1,tau2+1,tau3+1,:)=fft(r(tau1+1,tau2+1,tau3+1,:))/T;
end
end
end
cyclostationary_toolbox/cyclic_autocorrelation.m 0100644 0031654 0031654 00000001176 06436272744 0022750 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_autocorrelation(x,alpha,max_tau)
%
% CYCLIC_AUTOCORRELATION
%
% calculates the cyclic autocorrelation for signal
% x at frequency alpha
%
% R(k*alpha,tau)=E{x(t-tau/2)x(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
%
% USAGE
% R=cyclic_autocorrelation(x,alpha,max_tau)
%
% calculate cross correlation up to max_tau time lags
% File: cyclic_autocorrelation.m
% Last Revised: 24/11/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
R=cyclic_cross_correlation(x,x,alpha,max_tau); cyclostationary_toolbox/cyclic_autocorrelation_fast.m 0100644 0031654 0031654 00000002131 06436273635 0023755 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_autocorrelation_fast(x,T,max_tau)
%
% CYCLIC_AUTOCORRELATION_FAST
%
% calculates the cyclic autocorrelation for signal x
% at frequency alpha=1/T using a fast
% approximation based on the synchronous average of the
% time varying autocorrelation. Fundamental signal
% period can be defined as a single period or as a sequence
% of once per period pulse times.
%
% R(k*alpha,tau)=E{x(t-tau/2)x(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
%
% R=cyclic_autocorrelation_fast(x,T,max_tau)
%
% calculate autocorrelation up to max_tau time lags
%
% if T is a scalar, then T defined the fundamental
% cyclic period
%
% if T is a vector, then T defines a series of once
% per revolution impulses
% File: cyclic_autocorrelation_fast.m
% Last Revised: 24/11/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
R=cyclic_cross_correlation_fast(x,x,T,max_tau); cyclostationary_toolbox/cyclic_autocovariance.m 0100644 0031654 0031654 00000001256 06436273252 0022533 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_autocovariance(x,alpha,max_tau)
%
% CYCLIC_AUTOCOVARIANCE
%
% calculates the cyclic autocovariance for signal
% x at frequency alpha
%
% R(k*alpha,tau)=E{x(t-tau/2)x(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
% Signal x has its periodic average removed
%
% USAGE
% R=cyclic_autocovariance(x,alpha,max_tau)
%
% calculate autocovariance up to max_tau time lags
% File: cyclic_autocovariance.m
% Last Revised: 24/11/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
R=cyclic_cross_covariance(x,x,alpha,max_tau); cyclostationary_toolbox/cyclic_autocovariance_fast.m 0100644 0031654 0031654 00000002277 06506437754 0023563 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_autocovariance_fast(x,T,max_tau)
%
% CYCLIC_AUTOCOVARIANCE_FAST
%
% calculates the cyclic autocovariance for signal x
% at frequency alpha=1/T using a fast
% approximation based on the synchronous average of the
% time varying autocorrelation. Fundamental signal
% period can be defined as a single period or as a sequence
% of once per period pulse times.
%
% R(k*alpha,tau)=E{ x(t-tau/2) y(t+tau/2)
% exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
% Signals x has its periodic average removed
%
% USAGE
% R=cyclic_autocovariance_fast(x,T,max_tau)
%
% calculate autocorrelation up to max_tau time lags
%
% if T is a scalar, then T defined the fundamental
% cyclic period
%
% if T is a vector, then T defines a series of once
% per revolution impulses
% File: cyclic_autocovariance_fast.m
% Last Revised: 24/11/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
R=cyclic_cross_covariance_fast(x,x,T,max_tau);
cyclostationary_toolbox/cyclic_correlation_coefficient.m 0100644 0031654 0031654 00000002027 06436535324 0024405 0 ustar 00acmc acmc 0000304 0001726 function rho=cyclic_correlation_coefficient(S,tol)
%
% CYCLIC_CORRELATION_COEFFICIENT
%
% calculates the cyclic correlation coefficient for signal
% from its spectral correlation density
%
% USAGE
% rho=cyclic_correlation_coefficient(S,tol)
%
% Optional parameter tol can be used to set all
% rho(x,y) to zero if S(x,y) is less than tol% of max(S)
% File: cyclic_correlation_coefficient.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
[rows,cols]=size(S);
rho=zeros(rows,cols);
S=abs(S);
mxs=max(max(S));
scale_factor=2*cols/rows;
for f=1:rows
for a=1:cols
fplus=f+floor(a/scale_factor);
fminus=f-floor(a/scale_factor);
if fplus>rows
fplus=fplus-rows;
end
if fminus<1
fminus=fminus+rows;
end
rho(f,a)=S(f,a)/sqrt(S(fplus,1)*S(fminus,1));
end
end
if nargin==2
[i,j]=find(S<tol*mxs/100);
for k=1:length(i)
rho(i(k),j(k))=0;
end
end
cyclostationary_toolbox/cyclic_cross_correlation.m 0100644 0031654 0031654 00000002545 06517616215 0023264 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_cross_correlation(x,y,alpha,max_tau)
%
% CYCLIC_CROSS_CORRELATION
%
% calculates the cyclic cross correlation between
% two signals x,y at frequency alpha
%
% R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... 2*pi/alpha-1
%
% USAGE
% R=cyclic_cross_correlation(x,y,alpha,max_tau)
%
% calculate cross correlation up to max_tau time lags
% File: cyclic_cross_correlation.m
% Last Revised: 23/4/98
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=4
error('Incorrect number of arguments for function cyclic_cross_correlation');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2 pi in function cyclic_cross_correlation');
end
T=ceil(2*pi/alpha)-1;
lx=length(x);
t=0:lx-1;
R=zeros(max_tau*2+1,T+1);
% Compute even time shift segments
for tau=-max_tau:2:max_tau
for k=0:T
R(tau+1+max_tau,k+1)=mean(x(1:lx-max_tau-tau).*y(max_tau+tau+1:lx) ...
.*exp(-j*k*alpha*t(1+(max_tau+tau)/2:lx-(max_tau+tau)/2)));
end
end
% Compute odd time shift segments
t=t+0.5;
for tau=-max_tau+1:2:max_tau
for k=0:T
R(tau+1+max_tau,k+1)=mean(x(1:lx-tau-max_tau).*y(max_tau+tau+1:lx) ...
.*exp(-j*k*alpha*t(1+(max_tau+tau-1)/2:lx-(max_tau+tau+1)/2)));
end
end
cyclostationary_toolbox/cyclic_cross_correlation_fast.m 0100644 0031654 0031654 00000005317 06536756011 0024301 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_cross_correlation_fast(x,y,T,max_tau)
%
% CYCLIC_CROSS_CORRELATION_FAST
%
% calculates the cyclic cross correlation between
% two signals x,y at frequency alpha=1/T using a fast
% approximation based on the synchronous average of the
% time varying autocorrelation. Fundamental signal
% period can be defined as a single period or as a sequence
% of once per period pulse times.
%
% R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
%
% R=cyclic_cross_correlation_fast(x,y,T,max_tau)
%
% calculate cross correlation up to max_tau time lags
%
% if T is a scalar, then T defined the fundamental
% cyclic period
%
% if T is a vector, then T defines a series of once
% per revolution impulses
% File: cyclic_cross_correlation_fast.m
% Last Revised: 23/4/98
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=4
error('Incorrect number of arguments for function cyclic_cross_correlation_fast');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function cyclic_cross_correlation_fast');
end
% Ensure that vectors are the right sizes and shapes
if length(x)>length(y)
x=x(1:length(y));
end
[rows,cols]=size(x);
if rows>cols
x=x';
end
[rows,cols]=size(y);
if rows>cols
y=y';
end
% Remove excess samples due to non-integer sampling
if length(T)==1
% remove jitter samples if non-integer T
if T~=floor(T)
cp=1;np=1;
while cp+T<length(x)
cp=cp+floor(T);
np=np+T;
if (np-cp)>1
x=[x(1:cp-1) x(cp+1:length(x))];
y=[y(1:cp-1) y(cp+1:length(y))];
np=np-1;
end
end
end
n=floor((length(x)-2*max_tau-1)/T);
else
% extract time series correlated with periodic pulses
rot_positions=T;
T=floor(median(diff(rot_positions)));
nx=[];
ny=[];
n=length(rot_positions)-2;
for k=1:n;
cp=rot_positions(k);
nx=[nx x(cp:cp+T-1)];
ny=[ny y(cp:cp+T-1)];
end
nx=[nx x(rot_positions(n+1):rot_positions(n+1)+2*max_tau+1)];
x=nx;
ny=[ny y(rot_positions(n+1):rot_positions(n+1)+2*max_tau+1)];
y=ny;
end
% Compute time varying cross correlation
r=zeros(2*max_tau+1,floor(T));
temp=zeros(floor(T),n);
t=(1:floor(T)*n);
for k=-max_tau:max_tau
temp(:)=x(t+max_tau).*y(t+k+max_tau);
r(k+1+max_tau,:)=mean(temp');
end
% Take DFT of time varying correlation with appropriate phase change
% to compensate for time shift
R=zeros(2*max_tau+1,floor(T));
for k=-max_tau:max_tau
R(k+1+max_tau,:)=exp(-i*pi*((0:floor(T)-1)/T)*k).*fft(r(k+1+max_tau,:))/T;
end
cyclostationary_toolbox/cyclic_cross_covariance.m 0100644 0031654 0031654 00000003204 06517616452 0023051 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_cross_covariance(x,y,alpha,max_tau)
%
% CYCLIC_CROSS_COVARIANCE
%
% calculates the cyclic cross covariance between
% two signals x,y at frequency alpha
%
% R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
% With signals x and y having their periodic average removed
%
% USAGE
% R=cyclic_cross_covariance(x,y,alpha,max_tau)
%
% calculate cross covariance up to max_tau time lags
% File: cyclic_cross_covariance.m
% Last Revised: 23/4/98
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=4
error('Incorrect number of arguments for function cyclic_cross_covariance');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2 pi in function cyclic_cross_covariance');
end
% Remove cyclic mean from signals
cmx=cyclic_mean(x,alpha);
cmy=cyclic_mean(y,alpha);
lx=length(x);
t=0:lx-1;
T=ceil(2*pi/alpha)-1;
for k=1:lx
x(k)=x(k)-1/(2*pi)*sum(cmx.*exp(j*alpha*(0:T)*(k-1)));
y(k)=y(k)-1/(2*pi)*sum(cmy.*exp(j*alpha*(0:T)*(k-1)));
end
R=zeros(2*max_tau,T+1);
% Compute even time shift segments
for tau=-max_tau:2:max_tau
for k=0:T
R(tau+1+max_tau,k+1)=mean(x(1:lx-tau-max_tau).*y(tau+1+max_tau:lx) ...
.*exp(j*k*alpha*t(1+(tau+max_tau)/2:lx-(tau+max_tau)/2)));
end
end
% Compute odd time shift segments
t=t+0.5;
for tau=-max_tau+1:2:max_tau
for k=0:T
R(tau+1+max_tau,k+1)=mean(x(1:lx-tau-max_tau).*y(tau+1+max_tau:lx) ...
.*exp(j*k*alpha*t(1+(tau+max_tau-1)/2:lx-(tau+max_tau+1)/2)));
end
end
ormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=4
error('Incorrect number of arguments for function cyclic_cross_covariance');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2 pi in function cyclic_cross_covariance');
end
% Remove cyclic mean from signals
cmx=cyclic_mean(x,alpha);
cmy=cyclic_mean(y,alpha);
lx=length(x);
t=0:lx-1;
T=cyclostationary_toolbox/cyclic_cross_covariance_fast.m 0100644 0031654 0031654 00000006022 06544166221 0024061 0 ustar 00acmc acmc 0000304 0001726 function R=cyclic_cross_covariance_fast(x,y,T,max_tau)
%
% CYCLIC_CROSS_COVARIANCE_FAST
%
% calculates the cyclic cross covariance between
% two signals x,y at frequency alpha=1/T using a fast
% approximation based on the synchronous average of the
% time varying autocorrelation. Fundamental signal
% period can be defined as a single period or as a sequence
% of once per period pulse times.
%
% R(k*alpha,tau)=E{ x(t-tau/2) y(t+tau/2)
% exp(-jk(alpha)t)}
% for k=0 ... 1/alpha
% With signals x and y having their periodic average removed
%
% USAGE
% R=cyclic_cross_covariance_fast(x,y,T,max_tau)
%
% calculate cross correlation up to max_tau time lags
%
% if T is a scalar, then T defined the fundamental
% cyclic period
%
% if T is a vector, then T defines a series of once
% per revolution impulses
% File: cyclic_cross_covariance_fast.m
% Last Revised: 23/4/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=4
error('Incorrect number of arguments for function cyclic_cross_covariance_fast');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function cyclic_cross_covariance_fast');
end
% Ensure that vectors are the right sizes and shapes
if length(x)>length(y)
x=x(1:length(y));
end
[rows,cols]=size(x);
if rows>cols
x=x';
end
[rows,cols]=size(y);
if rows>cols
y=y';
end
% Calculate synchronous averages from signals
mx=synchronous_average(x,T);
my=synchronous_average(y,T);
% Remove excess samples due to non-integer sampling
% and renove cyclic mean from signal
if length(T)==1
cp=1;np=1;
while cp+T<length(x)
x(cp:cp+floor(T)-1)=x(cp:cp+floor(T)-1)-mx;
y(cp:cp+floor(T)-1)=y(cp:cp+floor(T)-1)-my;
cp=cp+floor(T);
np=np+T;
if (np-cp)>1
x=[x(1:cp-1) x(cp+1:length(x))];
y=[y(1:cp-1) y(cp+1:length(y))];
np=np-1;
end
end
n=floor((length(x)-2*max_tau-1)/T);
else
% extract time series correlated with periodic pulses
rot_positions=T;
T=floor(median(diff(rot_positions)));
nx=[];
ny=[];
n=length(rot_positions)-2;
for k=1:n;
cp=rot_positions(k);
nx=[nx; x(cp:cp+T-1)-mx];
ny=[ny; y(cp:cp+T-1)-my];
end
nx=[nx;x(rot_positions(n+1):rot_positions(n+1)+max_tau+1)-mx(1:max_tau+1)];
x=nx;
ny=[ny;y(rot_positions(n+1):rot_positions(n+1)+max_tau+1)-my(1:max_tau+1)];
y=ny;
end
% Compute time varying cross correlation
r=zeros(2*max_tau+1,floor(T));
temp=zeros(floor(T),n);
t=(1:floor(T)*n);
for k=-max_tau:max_tau
temp(:)=x(t+max_tau).*y(t+k+max_tau);
r(k+1+max_tau,:)=mean(temp');
end
% Take DFT of time varying correlation with appropriate phase change
% to compensate for time shift
R=zeros(2*max_tau+1,floor(T));
for k=-max_tau:max_tau
R(k+1+max_tau,:)=exp(-i*pi*((0:floor(T)-1)/T)*k).*fft(r(k+1+max_tau,:))/T;
end
cyclostationary_toolbox/cyclic_cross_periodogram.m 0100644 0031654 0031654 00000002633 06510420350 0023233 0 ustar 00acmc acmc 0000304 0001726 function S=cyclic_cross_periodogram(x,y,N,a,g,L)
%
% CYCLIC_CROSS_PERIODOGRAM
%
% calculate the spectral correlation density function
% of signals x and y using the cyclic cross periodogram
% algorithm
%
% Input parameters are:
% x,y signals
% N length of time window used for estimating frequency
% segments
% a window used for smoothing segments
% g window for smoothing correlation
% L decimation factor
%
% USAGE S=cyclic_cross_periodogram(x,y,N,a,g,L)
%
% e.g s=cyclic_cross_periodogram(s1,s1,64,'hamming','hamming',1)
if ~exist('L')
L=1;
end
lx=length(x);
if (length(y)~=lx)
error('Time series x and y must be same length')
end
n=0:floor((lx-N)/L);
ln=length(n);
a=feval(a,N)';
g=feval(g,ln)';
g=g/sum(g);
a=a/sum(a);
X=zeros(2*N+1,ln);
Y=zeros(2*N+1,ln);
Ts=1/N;
for f=-N:N
xf=x.*exp(-j*2*pi*f*(0:lx-1)*Ts);
yf=y.*exp(-j*2*pi*f*(0:lx-1)*Ts);
for i=1:ln
n_r=n(i)*L+(1:N);
X(f+N+1,i)=sum(a.*xf(n_r));
Y(f+N+1,i)=conj(sum(a.*yf(n_r)));
end
end
S=zeros(N+1,floor(N/2)+1);
for alpha=-floor(N/4):floor(N/4)
for f=-floor(N/2):floor(N/2)
f1=f+alpha;
f2=f-alpha;
if ((abs(f1)<N/2)&(abs(f2)<N/2))
S(f+floor(N/2)+1,floor(N/4)+alpha+1)=sum(g.*X(f1+N+1,:).*Y(f2+N+1,:));
end
end
end
cyclostationary_toolbox/cyclic_cumulants_fast.m 0100644 0031654 0031654 00000003075 06436526626 0022567 0 ustar 00acmc acmc 0000304 0001726 function [m,v,s,k]=cyclic_cumulants_fast(x,T)
%
% CYCLIC_CUMULANTS_FAST
% calculate synchronously averaged zero lag cumulants
% up to fourth order
%
% USGAE
% [m,v,s,k]=cyclic_cumulants_fast(x,T)
%
% where x is a cyclostationary time series
% and T is the fundamental period of interest
% or a vector of period start positions
%
% output parameters are:
% m - synchronous average ( cyclic mean )
% v - cyclic variance
% s - cyclic skewness
% k - cyclic kurosis
% File: cyclic_cumulants_fast.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
%remove jitter samples for non-integer T
if length(T)==1
cp=1;np=1;
while cp+T<length(x)
cp=cp+floor(T);
np=np+T;
if (np-cp)>1
x=[x(1:cp-1);x(cp+1:length(x))];
np=np-1;
end
end
n=floor((length(x)-1)/T);
nts=x;
else
rot_positions=T;
T=floor(median(diff(rot_positions)));
nts=[];
n=length(rot_positions)-2;
for k=1:n;
cp=rot_positions(k);
nts=[nts; x(cp:cp+T-1)];
end
nts=[nts;x(rot_positions(n+1):rot_positions(n+1)+1)];
end
temp=zeros(floor(T),n);
z=(1:floor(T)*n);
% calculate cyclic mean - time domain average
temp(:)=nts(z);
m1=mean(temp');
m2=mean(temp'.^2);
m3=mean(temp'.^3);
m4=mean(temp'.^4);
m=m1;
v=m2-m1.^2;
s=(m3-3*m2.*m1+2*m1.^3)./m2.^(1.5);
k=(m4-3*m2.^2+4*m3.*m1+12*m2.*m1.^2-6*m1.^4)./m2.^2;
cyclostationary_toolbox/cyclic_mean.m 0100644 0031654 0031654 00000001253 06436565332 0020450 0 ustar 00acmc acmc 0000304 0001726 function m=cyclic_mean(x,alpha)
%
% CYCLIC_MEAN Calculates the cyclic-mean of a signal
% m(k*alpha)=E{x(t)exp(-j(k*alpha)t)}
% for k = 0 ... 1/alpha
%
% USAGE
% m=cyclic_mean(x,alpha)
%
% File: cyclic_mean.m
% Last Revised: 25/11/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=2
error('Incorrect number of arguments for function cyclic_mean');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2*pi in function cyclic_mean');
end
T=ceil(2*pi/alpha)-1;
t=0:length(x)-1;
m=zeros(1,T+1);
for k=0:T
m(k+1)=mean(x.*exp(-j*k*alpha*t));
end
cyclostationary_toolbox/degree_of_cyclostationarity.m 0100644 0031654 0031654 00000001237 06436536615 0023771 0 ustar 00acmc acmc 0000304 0001726 function d=degree_of_cyclostationarity(R)
%
% DEGREE_OF_CYCLOSTATIONARITY
% Compute the degrees of cyclostationarity of a signal
% from its cyclic autocorrelation
%
% USAGE
% d=degree_of_cyclostationarity(R)
%
% R is the cyclic-autocorrelation of a signal
%
% d is a vector containing all degrees of cyclostationarity
% for integer multiples of the frequency used to obtain R
% File: degree_of_cyclostationarity.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
d=sum(abs(R).^2)/sum(abs(R(:,1)).^2);
cyclostationary_toolbox/fft_accumulation_method.m 0100644 0031654 0031654 00000003547 07045266605 0023073 0 ustar 00acmc acmc 0000304 0001726 function S=fft_accumulation_method(x,y,N,a,g,L)
%
% FFT_ACCUMULATION_METHOD
% estimate the spectral correlation density using the FFT
% accumulation method
%
% Reference:
% Roberts R. S., Brown, W. A. and Loomis, H. H.
% "Computationally efficient algorithms for cyclic
% spectral analysis" IEEE Signal Processing Magazine
% 8(2) pp38-49 April 1991
%
% Input parameters are:
% x,y signals
% N length of time window used for estimating frequency
% segments (should be power of 2)
% a window used for smoothing segments
% g window for smoothing correlation
% L decimation factor
%
% USAGE S=fft_accumulation_method(x,y,N,a,g,L)
%
% e.g s=fft_accumulation_method(s1,s1,64,'hamming','hamming',1)
if ~exist('L')
L=1;
end
lx=length(x);
if (length(y)~=lx)
error('Time series x and y must be same length')
end
n=0:floor((lx-N)/L);
ln=length(n);
a=feval(a,N)';
g=feval(g,ln)';
g=g/sum(g);
a=a/sum(a);
X=zeros(N,ln);
Y=zeros(N,ln);
Ts=1/N;
for i=1:ln
n_r=n(i)*L+(1:N);
X(:,i)=fftshift(fft(a.*x(n_r)))';
Y(:,i)=fftshift(conj(fft(a.*y(n_r))))';
end
lnd2=floor(ln/2);
lnd4=floor(ln/4)+1;
ln3d4=lnd4+lnd2-2;
S=zeros(2*N-1,lnd2*(N+1));
for alpha=-N/2+1:N/2-1
for f=-N/2:N/2-1
f1=f+alpha;
f2=f-alpha;
if ((abs(f1)<N/2)&(abs(f2)<N/2))
f1=f1+N/2;
f2=f2+N/2;
fsh=fftshift(fft(g.*X(f1,:).*Y(f2,:)));
S(2*f+N,(alpha+N/2)*lnd2+(1:lnd2-1))=fsh(1,lnd4:ln3d4);
end
f1=f+alpha;
f2=f-alpha+1;
if ((abs(f1)<N/2)&(abs(f2)<N/2))
f1=f1+N/2;
f2=f2+N/2;
fsh=fftshift(fft(g.*X(f1,:).*Y(f2,:)));
S(2*f+N+1,(alpha+N/2)*lnd2+(1:lnd2-1)-lnd4+1)=fsh(1,lnd4:ln3d4);
end
end
end
cyclostationary_toolbox/frac_delay_cyclic_ac.m 0100644 0031654 0031654 00000004056 06523572022 0022257 0 ustar 00acmc acmc 0000304 0001726 function R=frac_delay_cyclic_ac(x,N,alpha,max_tau)
%
% FRAC_DELAY_CYCLIC_AC
%
% calculates the cyclic auto correlation for signal x
% at frequency alpha resampling 1/alpha to N samples
% max_tau is expressed at original sampling fequency
%
% R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}
% for k=0 ... N-1
%
% USAGE
% R=frac_delay_cyclic_ac(x,N,alpha,max_tau)
%
% calculate cyclic autocorrelation up to max_tau time lags
% File: frac_delay_cyclic_ac.m
% Last Revised: 5/5/98
% Created: 5/5/98
% Author: Andrew C. McCormick
% (C) University of Strathclyde
lx=length(x);
Rxx=zeros(lx-2*max_tau,2*max_tau+1);
for k=-max_tau:max_tau
Rxx(:,k+max_tau+1)=[x(max_tau+1:lx-max_tau).*x(k+(max_tau+1:lx-max_tau))]';
end
r=zeros(N,2*max_tau+1);
w=zeros(N,1);
T=1/alpha;
for k=1:lx-2*max_tau
cycle=floor(k/T);
position=(k-cycle*T)/T*N+1;
fp=floor(position);cp=ceil(position);
if (fp==cp)
if (fp==0)
w(N)=w(N)+1;
r(N,:)=r(N,:)+Rxx(k,:);
else
w(fp)=w(fp)+1;
r(fp,:)=r(fp,:)+Rxx(k,:);
end
else
if (fp==0)
w(N)=w(N)+position-fp;
r(N,:)=r(N,:)+Rxx(k,:)*(position-fp);
else
w(fp)=w(fp)+position-fp;
r(fp,:)=r(fp,:)+Rxx(k,:)*(position-fp);
end
if cp>N
w(1)=w(1)+cp-position;
r(1,:)=r(1,:)+Rxx(k,:)*(cp-position);
else
w(cp)=w(cp)+cp-position;
r(cp,:)=r(cp,:)+Rxx(k,:)*(cp-position);
end
end
end
% If insufficient samples to estimate r, fill in blanks with
% the neighbouring sample in r.
if (sum(w==0)>0)
warning('Insufficient samples for desired resolution');
for k=1:N
if (w(k)==0)
if k>1
r(k,:)=r(k-1,:);
end
else
r(k,:)=r(k,:)/w(k);
end
end
else
r=r./(w*ones(1,2*max_tau+1));
end
% Take DFT of time varying correlation with appropriate phase change
% to compensate for time shift
R=zeros(2*max_tau+1,N);
for k=-max_tau:max_tau
R(k+1+max_tau,:)=exp(-i*pi*((0:N-1)/N)*k).*fft(r(:,k+1+max_tau)');
end
cyclostationary_toolbox/get_impulses.m 0100644 0031654 0031654 00000001714 06436527631 0020704 0 ustar 00acmc acmc 0000304 0001726 function [p,m,md,s]=get_impulses(ts,threshold,max_std)
%
% GET_IMPULSES
% get impulses from a noisy series of pulses
%
% USAGE
% [p,m,md,s]=get_impulses(x,threshold,max_deviation)
%
% Detects pulses p as a change in the signal > threshold
% Changes closer than max_deviation percent of the
% median pulse difference are discarded.
%
% Provides mean(m), median (md) and standard deviation(s)
% of times between pulses
% File: get_impulses.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
d=diff(ts);
p=find(d>threshold);
md=median(diff(p));
x=find(diff(p)>md-(max_std*md/100));
p=p(x);
md=median(diff(p));
x=find(diff(p)<md+(max_std*md/100));
p=p(x);
dp=diff(p);
x=find(dp<md+(max_std*md/100));
md=median(dp(x));
m=mean(dp(x));
s=std(dp(x));
cyclostationary_toolbox/scd.m 0100644 0031654 0031654 00000001631 06517612640 0016746 0 ustar 00acmc acmc 0000304 0001726 function S=scd(R,window)
%
% SCD Compute the spectral correlation density of a signal
% from its cyclic autocorrelation
%
% USAGE
% S=scd(R,window)
%
% R is the cyclic-autocorrelation of a signal
% The optional parameter window can be used to select
% a window function, the default is a hamming window
% File: scd.m
% Last Revised: 23/4/98
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
if nargin==1
window='hamming';
end
[r,c]=size(R);
% convert even number of rows to odd number of rows
win=feval(window,r);
win=win*ones(1,c);
S=fft(R.*win);
if nargout==0
dispS=S(1:(r-1)/2,1:floor(c/2)+1);
x=(0:floor(c/2));
y=(1:(r-1)/2)-1;
x=x/c;y=y/r;
contour(x,y,abs(dispS))
title('Spectral Correlation Density');
xlabel('alpha * pi radians');
ylabel('f * pi radians');
end
cyclostationary_toolbox/strip_spectral_correlation.m 0100644 0031654 0031654 00000002530 07045266647 0023644 0 ustar 00acmc acmc 0000304 0001726 function S=strip_spectral_correlation(x,y,N,a,g)
%
% STRIP_SPECTRAL_CORRELATION
% estimate the spectral correlation density using the strip
% spectral correlation
% N.B. Produces a skewed SCDF in S whith axes
% f+alpha/2, alpha
%
% Reference:
% Roberts R. S., Brown, W. A. and Loomis, H. H.
% "Computationally efficient algorithms for cyclic
% spectral analysis" IEEE Signal Processing Magazine
% 8(2) pp38-49 April 1991
%
% Input parameters are:
% x,y signals
% N length of time window used for estimating frequency
% segments (should be power of 2)
% a window used for smoothing segments
% g window for smoothing correlation
%
% USAGE S=strip_spectral_correlation(x,y,N,a,g)
%
% e.g s=strip_spectral_correlation(s1,s1,64,'hamming','hamming')
lx=length(x);
ly=length(y);
ln=lx-N;
if (ln~=ly)
warning('Length y is not Length x+N, some data will not be used')
end
a=feval(a,N)';
g=feval(g,ln)';
g=g/sum(g);
a=a/sum(a);
X=zeros(N,ln);
for i=1:ln
n_r=(1:N)+i-1;
X(:,i)=fftshift(fft(a.*x(n_r)))';
end
l=min([ly ln]);
S=zeros(N,ln);
for f=1:N
S(2*f,:)=fftshift(fft(g.*X(f,1:l).*y(1:l)));
end
cyclostationary_toolbox/synchronous_average.m 0100644 0031654 0031654 00000002565 06544165456 0022300 0 ustar 00acmc acmc 0000304 0001726 function m=synchronous_average(x,T)
%
% SYNCHRONOUS_AVERAGE
% calculate the synchronous average of the signal x
% with period T
%
% USAGE
% m=synchronous_average(x,T)
%
% if T is a scalar, then T defined the fundamental
% cyclic period
%
% if T is a vector, then T defines a series of once
% per revolution impulses
% File: synchronous_average.m
% Last Revised: 24/11/97
% Created: 24/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
% Simple error checks
if nargin~=2
error('Incorrect number of arguments for function synchronous_average');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function synchronous average');
end
% Remove excess samples due to non-integer sampling
if length(T)==1
% remove jitter samples if non-integer T
if T~=floor(T)
cp=1;np=1;
while cp+T<length(x)
cp=cp+floor(T);
np=np+T;
if (np-cp)>1
x=[x(1:cp-1) x(cp+1:length(x))];
np=np-1;
end
end
end
n=floor((length(x)-1)/T);
else
% extract time series correlated with periodic pulses
rot_positions=T;
T=floor(median(diff(rot_positions)));
nx=[];
n=length(rot_positions)-2;
for k=1:n;
cp=rot_positions(k);
nx=[nx; x(cp:cp+T-1)];
end
x=nx;
end
temp=zeros(floor(T),n);
t=(1:floor(T)*n);
temp(:)=x(t);
m=mean(temp');
cyclostationary_toolbox/wvd.m 0100644 0031654 0031654 00000001160 06436522441 0016771 0 ustar 00acmc acmc 0000304 0001726 function W=wvd(S)
%
% WVD Compute the Wigner-Ville Time Frequency representaion of
% a signal from its spectral correlation density
%
% USAGE
% W=wvd(S);
%
% S is the scd of a signal
% File: wvd.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde
W=fft(S')';
[r,c]=size(W);
if nargout==0
dispW=(W(1:(r+1)/2,:));
x=(1:c)-1;
y=(1:(r+1)/2)-1;
y=y/r;
contour(x,y,abs(dispW))
title('Wigner-Ville Time-Frequency Distribution')
xlabel('Time/Samples')
ylabel('Frequency * pi radians')
end